module gpsbend_setup
  implicit none
  private
  public:: setup
        interface setup; module procedure setupbend; end interface

contains
subroutine setupbend(obsLL,odiagLL, &
        lunin,mype,awork,nele,nobs,toss_gps_sub,is,init_pass,last_pass,conv_diagsave)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    setupbend    compute rhs of oi for gps bending angle
!   prgmmr: cucurull, l.    org: JCSDA/NCEP           date: 2005-12-01
!
! abstract:  For gps bending angle observations, this routine
!              a) reads obs assigned to given mpi task (geographic region),
!              b) simulates obs from guess,
!              c) apply some quality control to obs,
!              d) load weight and innovation arrays used in minimization
!              e) collects statistics for runtime diagnostic output
!              f) writes additional diagnostic information to output file
!
! program history log:
!   2005-12-01  cucurull - original code
!   2005-12-09  derber   - remove psfcg and use ges_lnps instead
!   2005-12-21  treadon  - add super_gps, move some diagnostics statistics
!                         to genstat_gps
!   2006-01-04  treadon  - correct inconsistency when using qcfail
!   2006-02-02  treadon  - rename lnprsl as ges_lnprsl
!   2006-02-24  derber  - modify to take advantage of convinfo module
!   2006-02-24 cucurull - update QC parameters and compute representativeness error
!                       - fix bug when counting obs that fail gross qc check
!                       - fix bug when rejecting obs that fail sats QC
!   2006-04-14  middlecoff - changed IF test to avoid out-of-bounds-reference on DATA
!   2006-07-28  derber  - modify to use new inner loop obs data structure
!   2006-07-31  kleist  - change to use ges_ps instead of lnps
!   2006-09-20  cucurull - use geopotential heights at intermediate levels instead 
!                          of midpoint,levels,remove psges, generalize minimization terms
!                          to hybrid coordinate,penalize high level obs,new QC checks,
!                          remove obs above 30 km,add b_tkges, remove dtime,  
!                          improve obs pressure calculation for diagnostic purposes,
!                          increase the extended atmosphere from 6 to 10 levels
!   2006-10-20 cucurull - update QC statistical checks and representativeness error with 
!                         COSMIC data
!                       - add information to diagnostic file
!   2007-01-29 cucurull - remove observations below 6 km
!   2007-03-01 derber - add toss_gps_sub; simplify profile qc loop
!   2007-03-19 tremolet - binning of observations
!   2007-06-05 tremolet - add observation diagnostics structure
!   2007-04-18     su - add nchar and cdiagbuf to diagnostic file
!   2007-06-22 cucurull - generalize qc structure to enable regional GSI;
!                         reduce gpswork2;remove conv_diagsave from argument list;
!                         consistent data counts for qc checks;
!                         update diagnostic information to be consistent with other obs;
!                         modify diagnostic structure;
!                         fix bug for jprof in stats qc
!   2007-07-26 cucurull - update code to generalized vertical coordinate (3d pressure)
!   2007-09-21 cucurull - remove obs above 40km from qc checks
!   2008-04-14 treadon  - remove super_gps (moved to genstats_gps)
!   2008-05-23 safford  - rm unused vars and uses
!   2008-12-03 todling  - revisited Tremolet modifications in light of newer GSI
!                       - changed handle of tail%time
!   2009-08-19 guo      - changed for multi-pass setup with dtime_check().
!   2010-04-16 cucurull - substantial clean up, bugs fixes, and update according to ref    
!   2010-04-23 cucurull - simplify loops and define repe_gps here
!   2010-05-26 cucurull - modify ds
!   2010-06-11 cucurull - update Statistics QC
!   2010-05-24 guo      - remerged/reimplmented multi-pass setup in observer mode;
!   2010-08-09 lueken   - removed n_5km variable from code.
!   2010-08-10 treadon  - remove last check for gpshead allocate; clean up use statements, 
!                         replace (izero,ione) with (0,1), remove _i_kind suffix from integer 
!                         constants
!   2010-08-11 lcucurull - replace tpdpres with tpdpres(nobs) to fix bug in TL code
!   2010-08-18        hu - add tell to mpeu_util declaration
!   2010-10-25 cucurull  - add quality control options for C/NOFS satellite
!   2010-11-8  cucurull  - tune observational errors
!   2010-12-12 cucurull  - generalize number of layers above model top to nsig_ext
!   2011-01-05 cucurull  - add gpstop to reject anything above this value
!   2011-01-07 cucurull  - reject observation if outside new integration grid and compute
!                          the nsig_ext value the user needs to use
!   2011-01-13 lueken    - corrected init_pass and last_pass indentation
!   2011-01-18 cucurull - increase the size of mreal by one element to add gps_dtype information
!   2011-06-17 treadon  - remove call tell at end of routine
!   2011-08-16 cucurull - fix bug in statistics qc
!   2011-08-17 cucurull - add Oceansat-2, METOP-B GRAS, SAC-D, and M-T assimilation capabilities
!   2012-10-16 cucurull - add preliminary super-refraction QC, compute and include qrefqes in diagnostic file
!                         different meaning for toss_gps_sub, add Tandem-X assimilation capability
!   2013-01-26  parrish - change from grdcrd to grdcrd1, tintrp2a to tintrp2a1, tintrp2a11,
!                          tintrp3 to tintrp31 (to allow successful debug compile on WCOSS)
!   2013-10-19  todling - metguess now holds background
!   2014-04-10  todling - 4dvar fix: obs must be in current time bin
!   2014-12-30  derber - Modify for possibility of not using obsdiag
!   2015-10-01  guo   - full res obvsr: index to allow redistribution of obsdiags
!   2016-05-18  guo     - replaced ob_type with polymorphic obsNode through type casting
!   2016-06-24  guo     - fixed the default value of obsdiags(:,:)%tail%luse to luse(i)
!                       . removed (%dlat,%dlon) debris.
!   2016-11-29  shlyaeva - save linearized H(x) for EnKF
!   2017-02-09  guo     - Remove m_alloc, n_alloc.
!                       . Remove my_node with corrected typecast().
!   2019-08-21  Shao    - add COSMIC-2, metop-c and Paz 
!   2020-03-18  Shao    - update observation error for COSMIC-2
!   2020-04-13  Shao    - update the statistis QC for COSMIC-2
!   2020-05-21  Shao    - add comments to include commercial data ID information
!   2020-08-26  Shao/Bathmann - add Jacobian QC
!   2021-07-29  cucurull - revert gross error check to default values
!   2021-07-29  cucurull - fix forward operator issues identified with L127
!   2021-11-04  cucurull - turn off Jacbian QC
!   2021-11-05  cucurull - update QCs and optimize/improve forward operator; bug fixes
!   2022-01-28  cucurull - add Sentinel-6, PAZ
!   2022-04-06  collard  - reintroduce Jacbian QC as an option (default off)
!   2024-12-04  Li       - remove the QC check for rejecting MetOp data <8 km
!   2024-12-04  Li       - add GRACE-FO (803&804) data 
!   2024-12-04  Li       - add new obs error model by Chris Riedel
!   2025-09-08  Li       - add PlanetiQ YAM-8 (768) data
!
!   input argument list:
!     lunin    - unit from which to read observations
!     mype     - mpi task id
!     nele     - number of data elements per observation
!     nobs     - number of observations
!
!   output argument list:
!     awork    - array containing information for data counts and gross checks
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use mpeu_util, only: die,perr,tell,getindex
  use kinds, only: r_kind,i_kind
  use m_gpsStats, only: gps_allhead,gps_alltail
  use obsmod , only: nprof_gps,lobsdiag_allocated,&
      lobsdiagsave,nobskeep,&
      time_offset,lobsdiag_forenkf
  use qcmod, only: gps_jacqc
  use m_obsNode, only: obsNode
  use m_gpsNode , only: gpsNode
  use m_gpsNode , only: gpsNode_appendto
  use m_obsLList, only: obsLList
  use obsmod, only: luse_obsdiag
  use m_obsdiagNode, only: obs_diag
  use m_obsdiagNode, only: obs_diags
  use m_obsdiagNode, only: obsdiagLList_appendNode
  use m_obsdiagNode, only: obsdiagLList_nextNode
  use m_obsdiagNode, only: obsdiagNode_init
  use m_obsdiagNode, only: obsdiagNode_set
  use m_obsdiagNode, only: obsdiagNode_get
  use m_obsdiagNode, only: obsdiagNode_assert


  use gsi_4dvar, only: nobs_bins,hr_obsbin
  use guess_grids, only: ges_lnprsi,hrdifsig,geop_hgti,nfldsig
  use guess_grids, only: nsig_ext,gpstop,commgpstop,commgpserrinf
  use gridmod, only: nsig
  use gridmod, only: get_ij,latlon11
  use constants, only: fv,n_a,n_b,n_c,deg2rad,tiny_r_kind,r0_01,r18,r61,r63,r10000
  use constants, only: zero,half,one,two,eccentricity,semi_major_axis,&
      grav_equator,somigliana,flattening,grav_ratio,grav,rd,eps,three,four,five,&
      r100,r400
  use lagmod, only: setq, setq_TL
  use lagmod, only: slagdw, slagdw_TL
  use jfunc, only: jiter,miter,jiterstart
  use convinfo, only: cermin,cermax,cgross,cvar_b,cvar_pg,ictype
  use m_dtime, only: dtime_setup, dtime_check

  use m_gpsrhs, only: muse
  use m_gpsrhs, only: dbend_loc,xj
  use m_gpsrhs, only: n_t,n_q,n_p,nrefges
  use m_gpsrhs, only: rges,gp2gm,prsltmp_o,tges_o
  use m_gpsrhs, only: error,error_adjst
  use m_gpsrhs, only: ratio_errors
  use m_gpsrhs, only: rdiagbuf,cdiagbuf
  use m_gpsrhs, only: qcfail
  use m_gpsrhs, only: qcfail_jac
  use m_gpsrhs, only: qcfail_one,qcfail_two,qcfail_three,qcfail_five 
  use m_gpsrhs, only: qcfail_six,qcfail_seven 
  use m_gpsrhs, only: data_ier,data_igps,data_ihgt
  use m_gpsrhs, only: gpsrhs_alloc
  use m_gpsrhs, only: gpsrhs_dealloc
  use m_gpsrhs, only: gpsrhs_aliases
  use m_gpsrhs, only: gpsrhs_unaliases

  use state_vectors, only: levels, svars3d
  use gsi_bundlemod, only : gsi_bundlegetpointer
  use gsi_metguess_mod, only : gsi_metguess_get,gsi_metguess_bundle
  use sparsearr, only: sparr2, new, size, writearray
  implicit none

! Declare passed variables
  type(obsLList ),target,dimension(:),intent(in):: obsLL
  type(obs_diags),target,dimension(:),intent(in):: odiagLL
  integer(i_kind)                         ,intent(in   ) :: lunin,mype,nele,nobs
  real(r_kind),dimension(100+7*nsig)      ,intent(inout) :: awork
  real(r_kind),dimension(max(1,nprof_gps)),intent(inout) :: toss_gps_sub

  integer, intent(in):: is              ! index to GPSbend buffer variables
  logical, intent(in):: init_pass       ! flag the pass for the first background bin
  logical, intent(in):: last_pass       ! flag the pass for the last background bin
  logical, intent(in):: conv_diagsave   ! save diagnostics file

! Declare local parameters
  real(r_kind),parameter::  r240 = 240.0_r_kind
  real(r_kind),parameter:: six = 6.0_r_kind
  real(r_kind),parameter:: seven = 7.0_r_kind
  real(r_kind),parameter:: ten = 10.0_r_kind
  real(r_kind),parameter:: eight = 8.0_r_kind
  real(r_kind),parameter:: nine = 9.0_r_kind
  real(r_kind),parameter:: eleven = 11.0_r_kind
  real(r_kind),parameter:: r12=12.0_r_kind
  real(r_kind),parameter:: r20=20.0_r_kind
  real(r_kind),parameter:: r40=40.0_r_kind
  real(r_kind),parameter:: r80=80.0_r_kind
  real(r_kind),parameter:: r1em3 = 1.0e-3_r_kind
  real(r_kind),parameter:: r1em6 = 1.0e-6_r_kind
  character(len=*),parameter :: myname='setupbend'
  real(r_kind),parameter:: crit_grad = 157.0_r_kind
  real(r_kind),parameter:: r790000=790000.0_r_kind

! Declare local variables
  integer(i_kind):: grids_dim
  real(r_kind) cutoff,cutoff1,cutoff2,cutoff3,cutoff4,cutoff12,cutoff23,cutoff34
  real(r_kind) sin2,zsges,ds,ns
  real(r_kind),dimension(:),allocatable:: ddnj,grid_s,ref_rad_s

  real(r_kind) rsig,rsig_up,ddbend,tmean,qmean
  real(r_kind) termg,termr,termrg,hob,dbend,grad_mod
  real(r_kind) fact,pw,nrefges1,nrefges2,nrefges3,k4,delz
  real(r_kind) ratio,residual,obserror,obserrlm,cermaxuse,cerminuse,cgrossuse
  real(r_kind) errinv_input,errinv_adjst,errinv_final,err_final

  real(r_kind),dimension(nele,nobs):: data
  real(r_kind),dimension(nsig):: dbenddn,dbenddxi
  real(r_kind) pressure,hob_s,d_ref_rad,d_ref_rad_TL,hob_s_top
  real(r_kind),dimension(4) :: w4,dw4,dw4_TL
  
  integer(i_kind) ier,ilon,ilat,ihgt,igps,itime,ikx,iuse, &
                  iprof,ipctc,iroc,isatid,iptid,ilate,ilone,ioff,igeoid, &
                  ilsw,ilswflag
  integer(i_kind) i,j,k,kk,mreal,nreal,jj,ikxx,ibin
  integer(i_kind) mm1,nsig_up,ihob,istatus,nsigstart
  integer(i_kind) kprof,istat,k1,k2,nobs_out,top_layer_SR,bot_layer_SR,count_SR
  integer(i_kind),dimension(4) :: gps_ij
  integer(i_kind):: satellite_id,transmitter_id

  type(sparr2) :: dhx_dx
  integer(i_kind) :: iz, t_ind, q_ind, p_ind, nnz, nind

  real(r_kind),dimension(3,nsig+nsig_ext) :: q_w,q_w_tl
  real(r_kind),dimension(nsig) :: hges,irefges,zges,dhdt,dhdp
  real(r_kind),dimension(nsig+1) :: prsltmp
  real(r_kind),dimension(nsig,nsig)::dndp,dxidp
  real(r_kind),dimension(nsig,nsig)::dndt,dxidt,dndq,dxidq
  real(r_kind),dimension(nsig+nsig_ext) :: n_TL
  real(r_kind),dimension(0:nsig+nsig_ext+1) :: ref_rad,xi_TL
  real(r_kind),dimension(nsig+nsig_ext+20) :: ref_rad_out
  real(r_kind) :: dlat,dlon,rocprof,unprof,dtime,dpressure,trefges,qrefges
  real(r_kind) :: dbetan,dbetaxi,rdog,alt
  real(r_kind),dimension(nsig):: tges,qges,qges_o
  real(r_kind),dimension(nobs):: tpdpres

  logical,dimension(nobs):: luse
  integer(i_kind),dimension(nobs):: ioid ! initial (pre-distribution) obs ID
  logical proceed

  logical:: in_curbin, in_anybin, obs_check,qc_layer_SR, save_jacobian
  type(gpsNode),pointer:: my_head
  type(obs_diag),pointer:: my_diag
  type(obs_diags),pointer:: my_diagLL

  real(r_kind),allocatable,dimension(:,:,:  ) :: ges_z
  real(r_kind),allocatable,dimension(:,:,:,:) :: ges_tv
  real(r_kind),allocatable,dimension(:,:,:,:) :: ges_q

  type(obsLList),pointer,dimension(:):: gpshead
  logical:: commdat

  gpshead => obsLL(:)

  save_jacobian = conv_diagsave .and. jiter==jiterstart .and. lobsdiag_forenkf

!*******************************************************************************
! List of GPS RO satellites and corresponding BUFR id
!740 => COSMIC FM1
!741 => COSMIC FM2
!742 => COSMIC FM3
!743 => COSMIC FM4
!744 => COSMIC FM5
!745 => COSMIC FM6
!4   => MetOpA
!41  => Champ
!722 => GRACE A
!723 => GRACE B
!820 => SACC
!42  => TerraSAR-X
!43  => Tandem-X
!786 => C/NOFS
!421 => OCEANSAT-2
!3   => MetOpB
!440 => Megha-Tropiques
!821 => SACD
!44  => PAZ
!750-755 => COSMIC-2 Equatorial
!724-729 => COSMIC-2 Polar
!825 => KOMPSAT-5
!5   => MetOpC
!265 => GeoOptics CICERO OP1
!266 => GeoOptics CICERO OP2
!267 => PlanetiQ GNOMES-A
!268 => PlanetiQ GNOMES-B
!269 => Spire Lemur 3U CubeSat
!66 => Sentinel-6 
!803=> GRACE C
!804=> GRACE D

! Check to see if required guess fields are available
  call check_vars_(proceed)
  if(.not.proceed) return  ! not all vars available, simply return

! If require guess vars available, extract from bundle ...
  call init_vars_

! Read and reformat observations in work arrays.
  read(lunin)data,luse,ioid

  ier=1        ! index of obs error in data array (now 1/(original obs error))
  ilon=2       ! index of grid relative obs location (x)
  ilat=3       ! index of grid relative obs location (y)
  ihgt=4       ! index of obs vertical coordinate in data array - impact parameter
  igps=5       ! index of gps data (or residual) in data array
  itime=6      ! index of obs time in data array
  ikxx=7       ! index of observation type
  iprof=8      ! index of profile
  ipctc=9      ! index of percent confidence
  iroc=10      ! index of local radius of curvature (m)
  isatid=11    ! index of satellite identifier
  iptid=12     ! index of platform transmitter id number
  iuse=13      ! index of use parameter
  ilone=14     ! index of earth relative longitude (degrees)
  ilate=15     ! index of earth relative latitude (degrees)
  igeoid=16    ! index of geoid undulation (a value per profile, m) 
  ilsw=17      ! index of bending angle LSW 
  ilswflag=18  ! index of LSW flag 

! Intialize variables
  nsig_up=nsig+nsig_ext ! extend nsig_ext levels above interface level nsig
  rsig=real(nsig,r_kind)
  rdog=rd/grav
  rsig_up=real(nsig_up,r_kind)
  nobs_out=0
  hob_s_top=one
  mm1=mype+1
  ns=nsig/two
  nsigstart=nint(ns)
  ns=r80
  grids_dim=nint(ns)  ! grid points for integration of GPS bend
  ds=r10000
  allocate(ddnj(grids_dim),grid_s(grids_dim),ref_rad_s(grids_dim)) 

! Allocate arrays for output to diagnostic file
  mreal=25
  nreal=mreal
  if (lobsdiagsave) nreal=nreal+4*miter+1
  if (save_jacobian) then
    nnz = nsig * 3         ! number of non-zero elements in dH(x)/dx profile
    nind   = 3             ! number of dense subarrays 
    call new(dhx_dx, nnz, nind)
    nreal = nreal + size(dhx_dx)
    ! jacobian sparse array indices are the same for all obs and can be filled
    ! in once here:
    t_ind = getindex(svars3d, 'tv')
    if (t_ind < 0) then
      print *, 'Error: no variable tv in state vector. Exiting.'
      call stop2(1300)
    endif
    q_ind = getindex(svars3d, 'q')
    if (q_ind < 0) then
      print *, 'Error: no variable q in state vector. Exiting.'
      call stop2(1300)
    endif
    p_ind = getindex(svars3d, 'prse')
    if (p_ind < 0) then
      print *, 'Error: no variable prse in state vector. Exiting.'
      call stop2(1300)
    endif
    dhx_dx%st_ind(1)  = sum(levels(1:t_ind-1)) + 1
    dhx_dx%end_ind(1) = sum(levels(1:t_ind-1)) + nsig
    dhx_dx%st_ind(2)  = sum(levels(1:q_ind-1)) + 1
    dhx_dx%end_ind(2) = sum(levels(1:q_ind-1)) + nsig
    dhx_dx%st_ind(3)  = sum(levels(1:p_ind-1)) + 1
    dhx_dx%end_ind(3) = sum(levels(1:p_ind-1)) + nsig
  endif
  if(init_pass) call gpsrhs_alloc(is,'bend',nobs,nsig,nreal,grids_dim,nsig_ext)
  call gpsrhs_aliases(is)
  if(nreal/=size(rdiagbuf,1)) then
     call perr(myname,'unexpected dimension')
     call perr(myname,'nreal =',nreal)
     call perr(myname,'size(rdiagbuf,1) =',size(rdiagbuf,1))
     call die(myname)
  endif

  if(init_pass) then
     ! initialize explicitly workspace variables,
     ! although they have been initialized during the
     ! buffer_alloc() call.

     data_ier (:)=data(ier ,:)
     data_ihgt(:)=data(ihgt,:)
     data_igps(:)=data(igps,:)
     muse(:)=.false.

     !qc flags
     !qc = 1: beyond model top or bottom
     !qc = 2: gpstop (55 km or 45 km for commercial data)
     !qc = 3: gross check
     !qc = 4: super refraction
     !qc = 5: background bending angle > 0.05 rad
     !qc = 6: statistic (cutoff) QC
     !qc = 7: ddnj < 0
     !qcfail_jac: jacobian check (currently not used)
     qcfail=.false.
     qcfail_jac=zero


     qcfail_one=zero;qcfail_two=zero
     qcfail_three=zero;qcfail_five=zero
     qcfail_six=zero;qcfail_seven=zero
     toss_gps_sub=zero 
     dbend_loc=zero

  else ! (init_pass)

!    Restore those arrays saved from the previous pass

     data(ier ,:)=data_ier (:)
     data(ihgt,:)=data_ihgt(:)
     data(igps,:)=data_igps(:)
  endif ! (init_pass)

! define new equally spaced grid s
  do j=0,grids_dim-1
     grid_s(j+1)=j*ds
  enddo
  k4=n_c-n_a

! A loop over all obs.
  call dtime_setup()

  !$omp parallel do default(none), schedule(dynamic,1), &
  !$omp& firstprivate(itime,iuse,jiter,ilate,ilat,ilon,iroc,igeoid,ihgt,ikxx, &
  !$omp&   nsig,n_a,n_b,k4,ier,isatid,iptid,igps,iprof,gpstop,commgpstop, &
  !$omp&   deg2rad,mype,ilone,ilsw,ilswflag,nsig_up,grids_dim,rsig_up,ds, &
  !$omp&   eccentricity,tiny_r_kind), &
  !$omp& private(i,dtime,obs_check,in_curbin,in_anybin,sin2,dlat,dlon, &
  !$omp&   rocprof,unprof,ikx,prsltmp,tges,qges,hges,zsges,termg,termr, &
  !$omp&   termrg,qc_layer_SR,count_SR,top_layer_SR,bot_layer_SR,k,zges, &
  !$omp&   qmean,tmean,fact,pw,pressure,nrefges1,nrefges2,nrefges3, &
  !$omp&   irefges,ref_rad,qges_o,alt,grad_mod, &
  !$omp&   hob,satellite_id,transmitter_id, &
  !$omp&   kprof,dpressure,ihob,k1,k2,delz,trefges,qrefges, &
  !$omp&   commdat,d_ref_rad,q_w,ref_rad_s,hob_s,w4,dw4,ddnj,kk,ref_rad_out, &
  !$omp&   dbend,j,ddbend,cgrossuse,cermaxuse,cerminuse,obserror, &
  !$omp&   obserrlm,residual,ratio,cutoff,cutoff1,cutoff2,cutoff3,cutoff4, &
  !$omp&   cutoff12,cutoff23,cutoff34), &
  !$omp& shared(nobs,data,muse,tpdpres,ges_lnprsi,hrdifsig,nfldsig, &
  !$omp&   prsltmp_o,nrefges,grav,gp2gm,rges,tges_o,eps,luse,ictype, &
  !$omp&   n_q,n_p,n_t,fv,nsigstart,ges_tv,ges_q,geop_hgti,ges_z, &
  !$omp&   error,error_adjst,rsig,ratio_errors,qcfail_one,cdiagbuf, &
  !$omp&   rdiagbuf,qcfail,xj,dbend_loc,qcfail_three,cgross,time_offset,nsig_ext, &
  !$omp&   cermax,cermin,qcfail_seven,qcfail_five,qcfail_six,qcfail_two,grid_s), &
  !$omp& reduction(+:nobs_out,awork), reduction(max:toss_gps_sub,hob_s_top)
  loopoverobs1: do i=1,nobs ! loop over obs
     dtime=data(itime,i)
     obs_check=.false. 

     call dtime_check(dtime, in_curbin, in_anybin)
     if(.not.in_anybin) cycle   ! not interested, not even for ibin
     if(.not.in_curbin) cycle   ! skip rest of lookp

     muse(i)=nint(data(iuse,i)) <= jiter
     sin2  = sin(data(ilate,i)*deg2rad)**2
     dlat=data(ilat,i)
     dlon=data(ilon,i)
     rocprof=data(iroc,i)
     unprof=data(igeoid,i)
     tpdpres(i)=data(ihgt,i)
     ikx=nint(data(ikxx,i))

!    Interpolate log(pres),temperature,specific humidity, 
!    corrected geopotential heights and topography to obs location
     call tintrp2a1(ges_lnprsi,prsltmp,dlat,dlon,dtime,hrdifsig,nsig+1,mype,nfldsig)
     call tintrp2a1(    ges_tv,   tges,dlat,dlon,dtime,hrdifsig,nsig,mype,nfldsig)
     call tintrp2a1(     ges_q,   qges,dlat,dlon,dtime,hrdifsig,nsig,mype,nfldsig)
     call tintrp2a1( geop_hgti,   hges,dlat,dlon,dtime,hrdifsig,nsig+1,mype,nfldsig)
     call tintrp2a11(    ges_z,  zsges,dlat,dlon,dtime,hrdifsig,mype,nfldsig)

     prsltmp_o(1:nsig,i)=prsltmp(1:nsig) ! needed in minimization

! Compute refractivity index-radius product at interface
!
!    Convert geopotential height at layer midpoints to geometric height using
!    equations (17, 20, 23) in MJ Mahoney's note "A discussion of various
!    measures of altitude" (2001).  Available on the web at
!    http://mtp.jpl.nasa.gov/notes/altitude/altitude.html
!
!    termg  = equation 17
!    termr  = equation 21
!    termrg = first term in the denominator of equation 23
!    zges   = equation 23
 
     termg = grav_equator * &
             ((one+somigliana*sin2)/sqrt(one-eccentricity*eccentricity*sin2))
     termr = semi_major_axis / (one + flattening + grav_ratio - two*flattening*sin2)
     termrg = (termg/grav)*termr

     qc_layer_SR=.false.
     count_SR=0
     top_layer_SR=0
     bot_layer_SR=0

     !!dir$ ivdep
     do k=1,nsig 
        zges(k) = (termr*hges(k)) / (termrg-hges(k))  ! eq (23) at interface (topo corrected)
        gp2gm(k,i)= termr/(termrg-hges(k))+((termr*hges(k))/(termrg-hges(k))**2)
        rges(k,i) = zges(k) + zsges + unprof + rocprof   ! radius r_i

        if(k>1) then
           qmean=(qges(k)+qges(k-1))/two
           tmean=(tges(k)+tges(k-1))/two
        else
           qmean=qges(1)
           tmean=tges(1)
        endif
        tges_o(k,i)=tmean  ! needed in minimization
        qges_o(k)=qmean
        fact=(one+fv*qmean)
        pw=eps+qmean*(one-eps)
        pressure=ten*exp(prsltmp(k)) ! pressure of interface level in mb
        nrefges1=n_a*(pressure/tmean)*fact
        nrefges2=n_b*qmean*pressure*fact**2/(tmean**2*pw)
        nrefges3=k4*fact*qmean*pressure/(tmean*pw)
        nrefges(k,i)=nrefges1+nrefges2+nrefges3 ! refractivity N_i
        irefges(k)= one+(r1em6*nrefges(k,i))    ! index of refractivity n_i
        ref_rad(k)=irefges(k)*rges(k,i)         ! refractivity index-radius product x_i

!       Terms for the minimization
        n_q(k,i)= &
           n_a*(pressure/tmean)*fv+&
           (n_b/(tmean**2*pw))*pressure*fact*(fact+qmean*two*fv)-&
           (n_b/(tmean**2*pw**2))*pressure*qmean*fact**2*(one-eps)+&
           (k4*pressure/(tmean*pw))*(fv*qmean+fact)-&
           (k4/(tmean*pw**2))*fact*qmean*pressure*(one-eps)
        n_p(k,i)= &
           ten*((n_a/tmean)*fact+&
           (n_b/(tmean**2*pw))*qmean*fact**2+&
           (k4/(tmean*pw))*qmean*fact)
        n_t(k,i)= &
           -n_a*fact*(pressure/tmean**2)-n_b*qmean*fact**2*two*&
           (pressure/(tmean**3*pw))-&
           (k4/(tmean**2*pw))*fact*qmean*pressure

     end do
     alt=(tpdpres(i)-rocprof-zsges-unprof)*r1em3
     if (alt<= five) then 
        do k=nsigstart,1,-1 
!       check for model SR layer at obs location
           grad_mod=1000.0_r_kind*(nrefges(k+1,i)-nrefges(k,i))/(rges(k+1,i)-rges(k,i))
           if (abs(grad_mod)>= half*crit_grad) then  ! SR - likely, to be used in obs SR qc
              qc_layer_SR=.true.   !SR-likely layer detected
           endif
           if (abs(grad_mod) >= 0.75_r_kind*crit_grad) then  !relax to close-to-SR conditions
              count_SR=count_SR+1 ! layers of SR                                           
              if (count_SR > 1 ) then
                 bot_layer_SR=k
              else
                 top_layer_SR=k
                 bot_layer_SR=top_layer_SR
              endif
           endif
        end do 
     endif 

!    locate observation in model vertical grid
     hob=tpdpres(i)
     call grdcrd1(hob,ref_rad(1),nsig,1)
     data(ihgt,i)=hob

!    initialize error, error_adjst, and err_final
     error(i)=tiny_r_kind
     error_adjst(i)=tiny_r_kind

     if (hob<three .or. hob>rsig) then
        data(ier,i) = zero
        ratio_errors(i) = zero
        muse(i)=.false.
        qcfail_one(i) = one
     endif

!    Increment obs counter along with low and high obs counters
     if (luse(i))then
        awork(1)=awork(1)+one
        if(hob > rsig) awork(3)=awork(3)+one
        if(hob <  one) awork(2)=awork(2)+one
     endif

!    Save some diagnostic information
!    occultation identification
     !satellite_id         = data(isatid,i) ! receiver occ id
     !transmitter_id       = data(iptid,i)  ! transmitter occ id
     !!$omp critical
     !write(cdiagbuf(i),'(2(i4.4))') satellite_id,transmitter_id
     !!$omp end critical

     rdiagbuf(:,i)         = zero

     rdiagbuf(1,i)         = ictype(ikx)        ! observation type
     rdiagbuf(20,i)        = one                ! uses gps_ref (one = use of bending angle)
     rdiagbuf(2,i)         = data(iprof,i)      ! profile identifier
     rdiagbuf(3,i)         = data(ilate,i)      ! lat in degrees
     rdiagbuf(4,i)         = data(ilone,i)      ! lon in degrees
     rdiagbuf(7,i)         = tpdpres(i)-rocprof ! impact height in meters
!    rdiagbuf(7,i)         = tpdpres(i)         ! impact parameter in meters
     rdiagbuf(8,i)         = dtime-time_offset  ! obs time (hours relative to analysis time)
!    rdiagbuf(9,i)         = data(ipctc,i)      ! input bufr qc - index of per cent confidence
     rdiagbuf(9,i)         = zsges              ! model terrain (m)
     rdiagbuf(11,i)        = data(iuse,i)       ! data usage flag
     rdiagbuf(17,i)        = data(igps,i)       ! bending angle observation (radians)
     rdiagbuf(19,i)        = hob                ! model vertical grid (interface) if monotone grid
     rdiagbuf(22,i)        = 1.e+10_r_kind      ! spread (filled in by EnKF)
     rdiagbuf(23,i)        = (data(ilsw,i)/data(igps,i))*100._r_kind ! Fractional LSW 
     rdiagbuf(24,i)        = 0.0_r_kind         ! holder for STD4060 in hybrid error model 
     rdiagbuf(25,i)        = data(ilswflag,i)   ! LSW flag 

     if(ratio_errors(i) > tiny_r_kind)  then ! obs inside model grid

       if (alt <= five) then
          if (top_layer_SR >= 1) then ! SR exists for at least one layer. Check if obs is inside                                                                       
             if ((tpdpres(i)<ref_rad(top_layer_SR+1)) .and. &
                    (tpdpres(i)<ref_rad(bot_layer_SR))) then !obs below model SR/close-to-SR layer
                    qcfail(i)=.true.
             elseif (tpdpres(i)>ref_rad(top_layer_SR+5)) then ! obs above SR/close-to-SR layer
                    qcfail(i)=.false.
                    if(hob < top_layer_SR+1) then !location might be aliased to the lower section of the non-monotonicity 
                      hob = tpdpres(i)
                      call grdcrd1(hob,ref_rad(top_layer_SR+1),nsig-top_layer_SR-1,1) ! only non-monotonic section above SR layer 
                      data(ihgt,i) = hob+top_layer_SR
                      hob = hob+top_layer_SR
                      rdiagbuf(19,i) = hob
                    endif
             else ! obs inside model SR/shadow or close-to-SR layer                                                                                                     
                    qcfail(i)=.true.
             endif
          endif

!         check for SR in obs, will be updated in genstats. 
          if ( data(igps,i) >= 0.03_r_kind .and. qc_layer_SR) then
             kprof = data(iprof,i)
             toss_gps_sub(kprof) = max (toss_gps_sub(kprof),data(igps,i))
          endif
       endif 

       alt=(tpdpres(i)-rocprof)*r1em3
!      get pressure (in mb), temperature and moisture at obs location
       call tintrp31(ges_lnprsi(:,:,1:nsig,:),dpressure,dlat,dlon,hob,&
              dtime,hrdifsig,mype,nfldsig)
       ihob=hob
       k1=min(max(1,ihob),nsig)
       k2=max(1,min(ihob+1,nsig))
       delz=hob-real(k1,r_kind)
       delz=max(zero,min(delz,one))
       trefges=tges_o(k1,i)*(one-delz)+tges_o(k2,i)*delz
       qrefges=qges_o(k1)*(one-delz)+qges_o(k2)*delz !Lidia

       rdiagbuf( 6,i)  = ten*exp(dpressure) ! pressure at obs location (hPa) if monotone grid
       rdiagbuf(18,i)  = trefges ! temperature at obs location (Kelvin) if monotone grid
       rdiagbuf(21,i)  = qrefges ! specific humidity at obs location (kg/kg) if monotone grid
       commdat=.false.
       if ( (data(isatid,i)>=265 .and. data(isatid,i)<=269) .or. (data(isatid,i)==768) ) commdat=.true.
       if (.not. qcfail(i)) then ! not SR

         ratio_errors(i) = data(ier,i)
  
         error(i)=one/data(ier,i) ! one/original error
         data(ier,i)=one/data(ier,i) ! one/original error
         error_adjst(i)= ratio_errors(i)* data(ier,i) !one/adjusted error

!        Extending atmosphere above interface level nsig
         d_ref_rad=ref_rad(nsig)-ref_rad(nsig-1)
         do k=1,nsig_ext
            ref_rad(nsig+k)=ref_rad(nsig)+ k*d_ref_rad ! extended x_i
            nrefges(nsig+k,i)=nrefges(nsig+k-1,i)**2/nrefges(nsig+k-2,i) ! exended N_i
         end do

!        Lagrange coefficients
         ref_rad(0)=ref_rad(3)
         ref_rad(nsig_up+1)=ref_rad(nsig_up-2)
         do k=1,nsig_up
            call setq(q_w(:,k),ref_rad(k-1:k+1),3)
         enddo

         muse(i)=nint(data(iuse,i)) <= jiter 

!        Get refractivity index-radius and [d(ln(n))/dx] in new grid.
         intloop: do j=1,grids_dim
           ref_rad_s(j)=sqrt(grid_s(j)*grid_s(j)+tpdpres(i)*tpdpres(i)) !x_j
           xj(j,i)=ref_rad_s(j)
           hob_s=ref_rad_s(j)
           call grdcrd1(hob_s,ref_rad(1),nsig_up,1)
           if(top_layer_SR >=1) then
             if(hob_s < top_layer_SR+1) then !correct if wrong location
                hob_s = ref_rad_s(j)
                call grdcrd1(hob_s,ref_rad(top_layer_SR+1),nsig_up-top_layer_SR-1,1)
                hob_s = hob_s+top_layer_SR
             endif
           endif
           dbend_loc(j,i)=hob_s  !location of x_j with respect to extended x_i
 

           if (hob_s < rsig_up) then  !obs inside the new grid
              ihob=hob_s

!             Compute refractivity and derivative at target points 
!             using Lagrange interpolators
              call slagdw(ref_rad(ihob-1:ihob+2),ref_rad_s(j),&
                   q_w(:,ihob),q_w(:,ihob+1),&
                   w4,dw4,4)
              if(ihob==1) then
                 w4(4)=w4(4)+w4(1); w4(1:3)=w4(2:4);w4(4)=zero
                 dw4(4)=dw4(4)+dw4(1);dw4(1:3)=dw4(2:4);dw4(4)=zero
                 ihob=ihob+1
              endif
              if(ihob==nsig_up-1) then
                 w4(1)=w4(1)+w4(4); w4(2:4)=w4(1:3);w4(1)=zero
                 dw4(1)=dw4(1)+dw4(4); dw4(2:4)=dw4(1:3);dw4(1)=zero
                 ihob=ihob-1
              endif
              ddnj(j)=dot_product(dw4,nrefges(ihob-1:ihob+2,i))!derivative (dN/dx)_j                                                                      
              if(ddnj(j)>zero) then
                 qcfail(i)=.true.
                 qcfail_seven(i) = one
                 data(ier,i) = zero
                 ratio_errors(i) = zero
                 muse(i)=.false.
                 cycle loopoverobs1 ! reject observation                                                                                                  
              endif
           else
              obs_check=.true.
              if (obs_check) then ! only once per obs here
                 nobs_out=nobs_out+1
                 d_ref_rad=ref_rad(nsig)-ref_rad(nsig-1)
                 do kk=1,20
                    ref_rad_out(nsig_up+kk)=ref_rad(nsig_up)+ kk*d_ref_rad ! extended x_i
                 end do
                 do kk=1,nsig_up
                    ref_rad_out(kk)=ref_rad(kk)
                 enddo
              endif
              hob_s=ref_rad_s(j)
              call grdcrd1(hob_s,ref_rad_out,nsig_up+20,1)
              hob_s_top=max(hob_s,hob_s_top) 
           endif !obs in new grid
         end do intloop

         if (obs_check) then      ! reject observation
            qcfail(i)=.true.
            data(ier,i) = zero
            ratio_errors(i) = zero
            muse(i)=.false.
            cycle loopoverobs1 
         endif

!        bending angle (radians)
         dbend=ds*ddnj(1)/ref_rad_s(1)
         do j=2,grids_dim
            ddbend=ds*ddnj(j)/ref_rad_s(j)
            dbend=dbend+two*ddbend
         end do
         dbend=-r1em6*tpdpres(i)*dbend  

!        Accumulate diagnostic information
         rdiagbuf( 5,i)  = (data(igps,i)-dbend)/data(igps,i) ! incremental bending angle (x100 %)

         data(igps,i)=data(igps,i)-dbend !innovation vector

!        Remove very large simulated values
         if(dbend > 0.05_r_kind) then
           data(ier,i) = zero
           ratio_errors(i) = zero
           qcfail_five(i) = one
           qcfail(i)=.true.
           muse(i)=.false.
         endif

         if (alt <= gpstop) then ! go into qc checks
            if ((alt <= commgpstop) .or. (.not.commdat)) then 
               cgrossuse=cgross(ikx)
               cermaxuse=cermax(ikx)
               cerminuse=cermin(ikx) 
!              Gross error check
               obserror = one/max(ratio_errors(i)*data(ier,i),tiny_r_kind)
               obserrlm = max(cerminuse,min(cermaxuse,obserror))
               residual = abs(data(igps,i))
               ratio    = residual/obserrlm

               if (ratio > cgrossuse) then
                   if (luse(i)) then
                      awork(4) = awork(4)+one
                   endif
                   qcfail_three(i) = one
                   data(ier,i) = zero
                   ratio_errors(i) = zero
                   muse(i)=.false.
               else   
!                  Statistics QC check if obs passed gross error check
                   cutoff=zero
                   if (((data(isatid,i) > 749).and.(data(isatid,i) < 756)).or.commdat) then
                      cutoff1=(-4.725_r_kind+0.045_r_kind*alt+0.005_r_kind*alt**2)*one/two
                   else
                      cutoff1=(-4.725_r_kind+0.045_r_kind*alt+0.005_r_kind*alt**2)*two/three
                   end if
                   cutoff2=1.5_r_kind+one*cos(data(ilate,i)*deg2rad)
                   if(trefges<=r240) then
                     cutoff3=two
                   else
                      cutoff3=0.005_r_kind*trefges**2-2.3_r_kind*trefges+266_r_kind
                   endif
                   if (((data(isatid,i) > 749).and.(data(isatid,i) < 756)).or.commdat) then
                      cutoff3=cutoff3*one/two
                   else
                      cutoff3=cutoff3*two/three
                   end if
                   if (((data(isatid,i) > 749).and.(data(isatid,i) < 756)).or.commdat) then
                      cutoff4=(four+eight*cos(data(ilate,i)*deg2rad))*one/two
                   else
                      cutoff4=(four+eight*cos(data(ilate,i)*deg2rad))*two/three
                   end if
                   cutoff12=((36_r_kind-alt)/two)*cutoff2+&
                            ((alt-34_r_kind)/two)*cutoff1
                   cutoff23=((eleven-alt)/two)*cutoff3+&
                            ((alt-nine)/two)*cutoff2
                   cutoff34=((six-alt)/two)*cutoff4+&
                            ((alt-four)/two)*cutoff3
                   if(alt>36_r_kind) cutoff=cutoff1
                   if((alt<=36_r_kind).and.(alt>34_r_kind)) cutoff=cutoff12
                   if((alt<=34_r_kind).and.(alt>eleven)) cutoff=cutoff2
                   if((alt<=eleven).and.(alt>nine)) cutoff=cutoff23
                   if((alt<=nine).and.(alt>six)) cutoff=cutoff3
                   if((alt<=six).and.(alt>four)) cutoff=cutoff34
                   if(alt<=four) cutoff=cutoff4

                   if (((data(isatid,i) > 749).and.(data(isatid,i) < 756)).or.commdat) then
                      cutoff=two*cutoff*r0_01
                   else
                      cutoff=three*cutoff*r0_01
                   end if
 
                   if(abs(rdiagbuf(5,i)) > cutoff) then
                      qcfail(i)=.true.
                      qcfail_six(i) = one
                      data(ier,i) = zero
                      ratio_errors(i) = zero
                      muse(i) = .false.
                   end if
               end if !gross qc check
            end if ! commdat < commgpstop
         end if ! alt < gpstop
!        Remove obs above gpstop or commgpstop  
         if((alt > gpstop) .or. (commdat .and. (alt > commgpstop))) then
           data(ier,i) = zero
           ratio_errors(i) = zero
           qcfail_two(i)=one
           muse(i)=.false.
         endif

       end if ! obs above super-refraction and shadow layers
     end if ! obs inside the vertical grid

  end do loopoverobs1 ! end of loop over observations
  !$omp end parallel do
  !write(6,'("setupbend: Number of obs considered and accepted " 2I10)') nobs, count(mask=muse .neqv. .false.)

  if (nobs_out>=1) then
     write(6,*)'WARNING GPSRO:',nobs_out,'obs outside integration grid. Increase nsig_ext to',&
     int(hob_s_top)-nsig+1
  endif

! Loop over observation profiles. Compute penalty
! terms, and accumulate statistics.
  if(last_pass) then

     do i=1,nobs

        if (qcfail(i) .or. qcfail_five(i) > zero .or. &
            qcfail_six(i) > zero .or. qcfail_seven(i) > zero) then
           data(ier,i) = zero
           ratio_errors(i) = zero
           muse(i) = .false.

!          Counting obs tossed due to stats qc check
!          This calculation will be updated in genstats_gps due to toss_gps_sub

           if (luse(i))then
              if(data(ilate,i)> r20) then
                 awork(22) = awork(22)+one                !NH
              else if(data(ilate,i)< -r20)then
                 awork(23) = awork(23)+one                !SH
              else
                 awork(24) = awork(24)+one                !TR
              end if
           end if
        end if
     end do
  endif ! (last_pass)

! Loop to load arrays used in statistics output
  call dtime_setup()
  do i=1,nobs
     dtime=data(itime,i)
     call dtime_check(dtime,in_curbin,in_anybin)
     if(.not.in_anybin) cycle

     if(last_pass) then
        if (ratio_errors(i)*data(ier,i) <= tiny_r_kind) muse(i) = .false.
        ikx=nint(data(ikxx,i))

 
        ! flags for observations that failed qc checks
        ! zero = observation is good
 
        if(qcfail(i))                 rdiagbuf(10,i) = four !modified in genstats due to toss_gps_sub
        if(qcfail_one(i) == one)      rdiagbuf(10,i) = one
        if(qcfail_two(i) == one)      rdiagbuf(10,i) = two
        if(qcfail_three(i) == one)    rdiagbuf(10,i) = three
        if(qcfail_five(i) == one)     rdiagbuf(10,i) = five
        if(qcfail_six(i) == one)      rdiagbuf(10,i) = six
        if(qcfail_seven(i) == one)    rdiagbuf(10,i) = seven

        if(muse(i)) then                    ! modified in genstats_gps due to toss_gps_sub
           rdiagbuf(12,i) = one             ! minimization usage flag (1=use, -1=not used)
        else
           rdiagbuf(12,i) = -one
        endif

        if (ratio_errors(i)*data(ier,i)>tiny_r_kind) then
           err_final = ratio_errors(i)*data(ier,i)
        else
           err_final = tiny_r_kind
        endif

        errinv_input  = tiny_r_kind
        errinv_adjst  = tiny_r_kind
        errinv_final  = tiny_r_kind


        if (error(i)>tiny_r_kind)       errinv_input=error(i)
        if (error_adjst(i)>tiny_r_kind) errinv_adjst=error_adjst(i)
        if (err_final>tiny_r_kind)      errinv_final=err_final

        rdiagbuf(13,i) = zero ! nonlinear qc relative weight - will be defined in genstats_gps
        rdiagbuf(14,i) = errinv_input ! original inverse gps obs error (rad**-1)
        rdiagbuf(15,i) = errinv_adjst ! original + represent error inverse gps 
                                      ! obs error (rad**-1)
        rdiagbuf(16,i) = errinv_final ! final inverse observation error due to 
                                      ! superob factor (rad**-1) and qc
                                      ! modified in genstats_gps
     endif ! (last_pass)

!    Link observation to appropriate observation bin
     if (nobs_bins>1) then
        ibin = NINT( dtime/hr_obsbin ) + 1
     else
        ibin = 1
     endif
     IF (ibin<1.OR.ibin>nobs_bins) write(6,*)mype,'Error nobs_bins, ibin=',nobs_bins,ibin
     if(luse_obsdiag) my_diagLL => odiagLL(ibin)

!    Link obs to diagnostics structure
     if(luse_obsdiag)then
        my_diag => obsdiagLList_nextNode(my_diagLL, &
                create=.not.lobsdiag_allocated, & ! either make-a-new or move-to-next
                idv=is, &
                iob=ioid(i), &
                ich=1, &
                elat=data(ilate,i), &
                elon=data(ilone,i), &
                luse=luse(i), &
                miter=miter)
        if (.not.associated(my_diag)) call die(myname,'a null obsdiagLList_nextNode, create =',.not.lobsdiag_allocated)
     endif

     if(last_pass) then
        if (nobskeep>0.and.luse_obsdiag) call obsdiagNode_get(my_diag, jiter=nobskeep, muse=muse(i))

!       Save values needed for generate of statistics for all observations
        if(.not. associated(gps_allhead(ibin)%head))then
           gps_allhead(ibin)%n_alloc = 0
           allocate(gps_allhead(ibin)%head,stat=istat)
           if(istat /= 0)write(6,*)' failure to write gps_allhead '
           gps_alltail(ibin)%head => gps_allhead(ibin)%head
        else
           allocate(gps_alltail(ibin)%head%llpoint,stat=istat)
           if(istat /= 0)write(6,*)' failure to write gps_alltail%llpoint '
           gps_alltail(ibin)%head => gps_alltail(ibin)%head%llpoint
        end if
        gps_allhead(ibin)%n_alloc = gps_allhead(ibin)%n_alloc +1
        gps_alltail(ibin)%n_alloc = gps_allhead(ibin)%n_alloc

        gps_alltail(ibin)%head%idv=is
        gps_alltail(ibin)%head%iob=ioid(i)
        gps_alltail(ibin)%head%elat= data(ilate,i)
        gps_alltail(ibin)%head%elon= data(ilone,i)

        allocate(gps_alltail(ibin)%head%rdiag(nreal),stat=istatus)
        if (istatus/=0) write(6,*)'SETUPBEND:  allocate error for gps_alldiag, istatus=',istatus

        gps_alltail(ibin)%head%ratio_err= ratio_errors(i)
        gps_alltail(ibin)%head%obserr   = data(ier,i)
        gps_alltail(ibin)%head%dataerr  = data(ier,i)*data(igps,i)
        gps_alltail(ibin)%head%pg       = cvar_pg(ikx)
        gps_alltail(ibin)%head%b        = cvar_b(ikx)
        gps_alltail(ibin)%head%loc      = data(ihgt,i)
        gps_alltail(ibin)%head%kprof    = data(iprof,i)
        gps_alltail(ibin)%head%type     = data(ikxx,i)
        gps_alltail(ibin)%head%luse     = luse(i) ! logical
        gps_alltail(ibin)%head%muse     = muse(i) ! logical
        satellite_id         = data(isatid,i) ! receiver occ id
        transmitter_id       = data(iptid,i)  ! transmitter occ id
        write(cdiagbuf(i),'(2(i4.4))') satellite_id,transmitter_id
        gps_alltail(ibin)%head%cdiag    = cdiagbuf(i)

!       Fill obs diagnostics structure
        if (luse_obsdiag) then
           call obsdiagNode_set(my_diag,wgtjo=(data(ier,i)*ratio_errors(i))**2, &
                jiter=jiter,muse=muse(i),nldepart=data(igps,i) )
        endif

!       Load additional obs diagnostic structure
        ioff = mreal
        if (lobsdiagsave) then
           associate(odiag => my_diag )
              do jj=1,miter
                 ioff=ioff+1
                 if (odiag%muse(jj)) then
                    rdiagbuf(ioff,i) = one
                 else
                    rdiagbuf(ioff,i) = -one
                 endif
              enddo
              do jj=1,miter+1
                 ioff=ioff+1
                 rdiagbuf(ioff,i) = odiag%nldepart(jj)
              enddo
              do jj=1,miter
                 ioff=ioff+1
                 rdiagbuf(ioff,i) = odiag%tldepart(jj)
              enddo
              do jj=1,miter
                 ioff=ioff+1
                 rdiagbuf(ioff,i) = odiag%obssen(jj)
              enddo
           end associate  ! odiag
        endif

        ! if obs is not "acceptable" and jacobian is not computed, fill jacobian
        ! with zeros
        if (save_jacobian) then
           dhx_dx%val = 0._r_kind
           call writearray(dhx_dx, rdiagbuf(ioff+1:nreal,i))
        endif

! If obs is "acceptable", load array with obs info for use
! in inner loop minimization (int* and stp* routines)

        if (in_curbin .and. muse(i)) then

           allocate(my_head)
           call gpsNode_appendto(my_head,gpshead(ibin))

           my_head%idv = is
           my_head%iob = ioid(i)
           my_head%elat= data(ilate,i)
           my_head%elon= data(ilone,i)

           allocate(my_head%jac_t(nsig),my_head%jac_q(nsig), &
                    my_head%jac_p(nsig+1),my_head%ij(4,nsig),stat=istatus)
           if (istatus/=0) write(6,*)'SETUPBEND:  allocate error for gps_point, istatus=',istatus

           gps_alltail(ibin)%head%mmpoint  => my_head
 
!          Inizialize some variables
           dxidt=zero; dxidp=zero; dxidq=zero
           dndt=zero; dndq=zero; dndp=zero

!          Set (i,j) indices of guess gridpoint that bound obs location
           call get_ij(mm1,data(ilat,i),data(ilon,i),gps_ij,my_head%wij)
 
!$omp parallel do  schedule(dynamic,1) private(k,j,dhdt,dhdp)
           do k=1,nsig
 
              my_head%ij(1,k)=gps_ij(1)+(k-1)*latlon11
              my_head%ij(2,k)=gps_ij(2)+(k-1)*latlon11
              my_head%ij(3,k)=gps_ij(3)+(k-1)*latlon11
              my_head%ij(4,k)=gps_ij(4)+(k-1)*latlon11
 
              dhdp=zero; dhdt=zero       
              if(k > 1) then
                 do j=2,k
                    dhdt(j-1)= rdog*(prsltmp_o(j-1,i)-prsltmp_o(j,i))
                    dhdp(j)= dhdp(j)-rdog*(tges_o(j-1,i)/exp(prsltmp_o(j,i)))
                    dhdp(j-1)=dhdp(j-1)+rdog*(tges_o(j-1,i)/exp(prsltmp_o(j-1,i)))
                 end do
              end if
              if(k == 1)then
                 dndt(k,k)=dndt(k,k)+n_t(k,i)
                 dndq(k,k)=dndq(k,k)+n_q(k,i)
                 dndp(k,k)=dndp(k,k)+n_p(k,i)
              else
                 dndt(k,k)=dndt(k,k)+half*n_t(k,i)
                 dndt(k,k-1)=dndt(k,k-1)+half*n_t(k,i)
                 dndq(k,k)=dndq(k,k)+half*n_q(k,i)
                 dndq(k,k-1)=dndq(k,k-1)+half*n_q(k,i)
                 dndp(k,k)=n_p(k,i)
              end if
              irefges(k)=one+r1em6*nrefges(k,i)
              ref_rad(k)=irefges(k)*rges(k,i)
              do j=1,nsig
                 dxidt(k,j)=r1em6*rges(k,i)*dndt(k,j)+irefges(k)*gp2gm(k,i)*dhdt(j)
                 dxidq(k,j)=r1em6*rges(k,i)*dndq(k,j)
                 dxidp(k,j)=r1em6*rges(k,i)*dndp(k,j)+irefges(k)*gp2gm(k,i)*dhdp(j)
              end do
           end do
           d_ref_rad=ref_rad(nsig)-ref_rad(nsig-1)
           do k=1,nsig_ext
              ref_rad(nsig+k)=ref_rad(nsig) + k*d_ref_rad
           end do
           ref_rad(0)=ref_rad(3)
           ref_rad(nsig_up+1)=ref_rad(nsig_up-2)
!$omp parallel do  schedule(dynamic,1) private(kk,k,j,xi_TL,n_TL,q_w,q_w_tl,d_ref_rad_TL,ihob,dw4,dw4_TL,dbetaxi,dbetan)
           do kk=1,nsig
              xi_TL=zero
              xi_TL(kk)=one
              n_TL=zero
              n_TL(kk)=one
              q_w=zero
              q_w_TL=zero
              d_ref_rad_TL=xi_TL(nsig)-xi_TL(nsig-1)
              do k=1,nsig_ext
                 xi_TL(nsig+k)=xi_TL(nsig)+ k*d_ref_rad_TL
                 n_TL(nsig+k)=(two*nrefges(nsig+k-1,i)*n_TL(nsig+k-1)/nrefges(nsig+k-2,i))-&
                     (nrefges(nsig+k-1,i)**2/nrefges(nsig+k-2,i)**2)*n_TL(nsig+k-2)

              end do
              xi_TL(0)=xi_TL(3)
              xi_TL(nsig_up+1)=xi_TL(nsig_up-2)
              do k=1,nsig_up
                 call setq_TL(q_w(:,k),q_w_TL(:,k),ref_rad(k-1:k+1),xi_TL(k-1:k+1),3)
              enddo
              intloop2: do j=1,grids_dim
                 ihob=dbend_loc(j,i)

! Compute refractivity and derivative at target points using Lagrange interpolators
                 call slagdw_TL(ref_rad(ihob-1:ihob+2),xi_TL(ihob-1:ihob+2),xj(j,i),&
                            q_w(:,ihob),q_w_TL(:,ihob),&
                            q_w(:,ihob+1),q_w_TL(:,ihob+1),&
                            dw4,dw4_TL,4)
                 if(ihob==1) then
                    dw4(4)=dw4(4)+dw4(1);dw4(1:3)=dw4(2:4);dw4(4)=zero
                    dw4_TL(4)=dw4_TL(4)+dw4_TL(1);dw4_TL(1:3)=dw4_TL(2:4);dw4_TL(4)=zero
                    ihob=ihob+1
                 endif
                 if(ihob==nsig_up-1) then
                    dw4(1)=dw4(1)+dw4(4); dw4(2:4)=dw4(1:3);dw4(1)=zero
                    dw4_TL(1)=dw4_TL(1)+dw4_TL(4); dw4_TL(2:4)=dw4_TL(1:3);dw4_TL(1)=zero
                    ihob=ihob-1
                 endif
                 dbetaxi=(r1em6/xj(j,i))*dot_product(dw4_TL,nrefges(ihob-1:ihob+2,i))
                 dbetan =(r1em6/xj(j,i))*dot_product(dw4,n_TL(ihob-1:ihob+2))
                 if(j == 1)then
                    dbenddxi(kk)=dbetaxi
                    dbenddn(kk)=dbetan
                 else
                    dbenddxi(kk)=dbenddxi(kk)+two*dbetaxi
                    dbenddn(kk)=dbenddn(kk)+two*dbetan
                 end if
              end do intloop2
              dbenddxi(kk)=-dbenddxi(kk)*ds*tpdpres(i)
              dbenddn(kk)=-dbenddn(kk)*ds*tpdpres(i)
           end do
!$omp parallel do  schedule(dynamic,1) private(k,j)
           do k=1,nsig
              my_head%jac_t(k)=zero
              my_head%jac_q(k)=zero
              my_head%jac_p(k)=zero
              do j=1,nsig
                 my_head%jac_t(k)=my_head%jac_t(k)+dbenddxi(j)*dxidt(j,k)+ &
                                                   dbenddn(j) * dndt(j,k)
                 my_head%jac_q(k)=my_head%jac_q(k)+dbenddxi(j)*dxidq(j,k)+ &
                                                   dbenddn(j) * dndq(j,k)
                 my_head%jac_p(k)=my_head%jac_p(k)+dbenddxi(j)*dxidp(j,k)+ &
                                                   dbenddn(j) * dndp(j,k)
              end do

              if (gps_jacqc .and. ((abs(my_head%jac_t(k)) > 0.0016_r_kind).or.(abs(my_head%jac_q(k)) > 7.5_r_kind).or. &
                  (abs(my_head%jac_p(k)) > 0.004_r_kind))) qcfail_jac(i) = one

           end do 

           my_head%jac_p(nsig+1) = zero

           if (qcfail_jac(i) == one) then
              do k=1,nsig
                 my_head%jac_t(k) = zero
                 my_head%jac_q(k) = zero
                 my_head%jac_p(k) = zero
              end do
              ratio_errors(i) = zero
              data(ier,i) = zero
              rdiagbuf(12,i) = -one
              rdiagbuf(10,i) = six
           end if

           if (save_jacobian) then
              ! fill in the jacobian
              do iz = 1, nsig
                 dhx_dx%val(iz)        = my_head%jac_t(iz)
                 dhx_dx%val(iz+nsig)   = my_head%jac_q(iz)
                 dhx_dx%val(iz+2*nsig) = my_head%jac_p(iz)
              enddo

              call writearray(dhx_dx, rdiagbuf(ioff+1:nreal,i))
              ioff = ioff + size(dhx_dx)
           endif

           my_head%raterr2= ratio_errors(i)**2     
           my_head%res    = data(igps,i)
           my_head%err2   = data(ier,i)**2
           my_head%time   = data(itime,i)
           my_head%b      = cvar_b(ikx)
           my_head%pg     = cvar_pg(ikx)
           my_head%luse   = luse(i)

           if (luse_obsdiag) then
              call obsdiagNode_assert(my_diag, my_head%idv,my_head%iob,1,myname,'my_diag:my_head')
              my_head%diags => my_diag
           endif

           my_head => null()
        end if ! (in_curbin .and. muse=1)
        do j=1,nreal
           gps_alltail(ibin)%head%rdiag(j)= rdiagbuf(j,i)
        end do
        gps_alltail(ibin)%head%ratio_err= ratio_errors(i)
        gps_alltail(ibin)%head%obserr   = data(ier,i)
        gps_alltail(ibin)%head%dataerr  = data(ier,i)*data(igps,i)
        gps_alltail(ibin)%head%muse     = muse(i) ! logical
     endif ! (last_pass)
  end do ! i=1,nobs
  deallocate(ddnj,grid_s,ref_rad_s)
  ! Release memory of local guess arrays
  call final_vars_

  ! Save these arrays for later passes
  data_ier (:)=data(ier ,:)
  data_ihgt(:)=data(ihgt,:)
  data_igps(:)=data(igps,:)

  call gpsrhs_unaliases(is)
  if(last_pass) call gpsrhs_dealloc(is)

  return
  contains

  subroutine check_vars_ (proceed)
  logical,intent(inout) :: proceed
  integer(i_kind) ivar,istatus
! Check to see if required guess fields are available
  call gsi_metguess_get ('var::q', ivar, istatus )
  proceed=ivar>0
  call gsi_metguess_get ('var::z' , ivar, istatus )
  proceed=proceed.and.ivar>0
  call gsi_metguess_get ('var::tv', ivar, istatus )
  proceed=proceed.and.ivar>0
  end subroutine check_vars_ 

  subroutine init_vars_

  real(r_kind),dimension(:,:  ),pointer:: rank2=>NULL()
  real(r_kind),dimension(:,:,:),pointer:: rank3=>NULL()
  character(len=5) :: varname
  integer(i_kind) ifld,istatus

! If require guess vars available, extract from bundle ...
  if(size(gsi_metguess_bundle)==nfldsig) then
!    get z ...
     varname='z'
     call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus)
     if (istatus==0) then
         if(allocated(ges_z))then
            write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '
            call stop2(999)
         endif
         allocate(ges_z(size(rank2,1),size(rank2,2),nfldsig))
         ges_z(:,:,1)=rank2
         do ifld=2,nfldsig
            call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus)
            ges_z(:,:,ifld)=rank2
         enddo
     else
         write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
         call stop2(999)
     endif
!    get tv ...
     varname='tv'
     call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus)
     if (istatus==0) then
         if(allocated(ges_tv))then
            write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '
            call stop2(999)
         endif
         allocate(ges_tv(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig))
         ges_tv(:,:,:,1)=rank3
         do ifld=2,nfldsig
            call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus)
            ges_tv(:,:,:,ifld)=rank3
         enddo
     else
         write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
         call stop2(999)
     endif
!    get q ...
     varname='q'
     call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus)
     if (istatus==0) then
         if(allocated(ges_q))then
            write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '
            call stop2(999)
         endif
         allocate(ges_q(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig))
         ges_q(:,:,:,1)=rank3
         do ifld=2,nfldsig
            call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus)
            ges_q(:,:,:,ifld)=rank3
         enddo
     else
         write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
         call stop2(999)
     endif
  else
     write(6,*) trim(myname), ': inconsistent vector sizes (nfldsig,size(metguess_bundle) ',&
                 nfldsig,size(gsi_metguess_bundle)
     call stop2(999)
  endif
  end subroutine init_vars_

  subroutine final_vars_
    if(allocated(ges_q )) deallocate(ges_q )
    if(allocated(ges_tv)) deallocate(ges_tv)
    if(allocated(ges_z )) deallocate(ges_z )
  end subroutine final_vars_

end subroutine setupbend
end module gpsbend_setup
